library(here)
library(data.table)
library(qs)
library(lfe)
library(tidyverse)
library(ggplot2)
library(ggthemes)
library(texreg)

# get the fuctios
source(here("code", "functions.R"))

# Get the data
vox <- qread(here("data", "vox.qs")); setDT(vox)

vox$age <- log(vox$age)

# Traditions --------------------------------------------------------------

#Would you like a Switzerland that is modern or a Switzerland that protects its traditions?
#Higher is traditional
vox$protectTrad_std <- (vox$protectTrad - mean(vox$protectTrad, na.rm = T)) / sd(vox$protectTrad, na.rm = T)
mod_protectTrad <- felm(protectTrad_std ~ transBorder15 + freeBorder15 + income + female + age + education + employed + married | bfs19 + voxnr | 0 | bfs19, weights = vox %>% filter(BR == 1 & travelMUNmin <= 30) %>% pull(gewteil), data = vox %>% filter(BR == 1 & travelMUNmin <= 30))


# Calm and Order --------------------------------------------------------------
#Would you like a Switzerland where peace and order are little emphasized, or a Switzerland where peace and order are strongly emphasized?
#Higher is strongly emphasized

## Distribution of variables by region -----------------------------------------
vox$TR_three <- NULL
vox$TR_three[vox$border15 == 1 & vox$BR == 1] <- "0-15 Minutes"
vox$TR_three[vox$border30 == 1 & vox$BR == 1] <- "15-30 Minutes"

calm_order_trends <- vox %>% filter(complete.cases(calmOrder), complete.cases(gewteil), complete.cases(TR_three)) %>%
  group_by(year, TR_three) %>%
  summarise(weighted_mean_calmOrder = sum(calmOrder * gewteil) / sum(gewteil), .groups = 'drop') %>%
  ggplot(aes(x = year, y = weighted_mean_calmOrder, color = TR_three, fill = TR_three, shape = TR_three, linetype = TR_three)) + 
  stat_summary(fun.y = mean, geom = "point") + 
  geom_line() + 
  geom_vline(xintercept = c(1999, 2004), linetype = "dashed", color = "black") +
  ylab("Outcome: 1-6\n
       (1: Not Emphasized, 6: Strongly Emphasized)") +
  xlab("Year") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1994, 2016, 2), expand = expansion(mult = 0, add = 0.5)) +
  scale_color_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_fill_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_linetype_manual(values = c("solid", "dashed"), name = NULL, guide = "legend") +
  scale_shape_manual(values = c(16, 17), name = NULL, guide = "legend") +
  scale_y_continuous(expand = expand_scale(mult = 0, add = 0.025)) +
  annotate("text", x = 1996, y = 5, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 5, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2007, y = 5, label = "Free Movement", hjust = 0.5) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"),
        axis.title.y = element_text(vjust = 0)) +
  ggtitle("Trends in Views on Law and Order")

ggsave(plot = calm_order_trends, here("figures", "VOX_calmOrder_trends.pdf"), height = 5, width = 6)


## Regression -----------------------------------------
vox$calmOrder_std <- (vox$calmOrder - mean(vox$calmOrder, na.rm = T)) / sd(vox$calmOrder, na.rm = T)

mod_calmOrder <- felm(calmOrder_std ~ transBorder15 + freeBorder15 + income + female + age + education + employed + married | bfs19 + voxnr | 0 | bfs19, weights = vox %>% filter(BR == 1 & travelMUNmin <= 30) %>% pull(gewteil), data = vox %>% filter(BR == 1 & travelMUNmin <= 30))




# Environment --------------------------------------------------------------
#•	a91j: would you like Switzerland to prioritize environmental protection over economic prosperity or vice versa? 
# •	Values:
#   1)	Environmental protection
# 2)	2
# 3)	3
# 4)	4
# 5)	5
# 6)	Economic prosperity
# 7)	Don’t know / no answer

vox$enviroProtect_std <- (vox$enviroProtect - mean(vox$enviroProtect, na.rm = T)) / sd(vox$enviroProtect, na.rm = T)

mod_enviroProtect <- felm(enviroProtect_std ~ transBorder15 + freeBorder15 + income + female + age + education + employed + married | bfs19 + voxnr | 0 | bfs19, weights = vox %>% filter(BR == 1 & travelMUNmin <= 30) %>% pull(gewteil), data = vox %>% filter(BR == 1 & travelMUNmin <= 30))


# Plot Main Paper --------------------------------------------------------------------
df_plot <- coef.gg.prep(list(mod_protectTrad)) %>% mutate(outcome = "Protect traditions") %>% 
  bind_rows(coef.gg.prep(list(mod_calmOrder)) %>% mutate(outcome = "Emphasize law & order"),
            coef.gg.prep(list(mod_enviroProtect)) %>% mutate(outcome = "Prosperity over environment")
            )

plot_culture_security <- ggplot(df_plot, aes(x = modelcoef, y = outcome, xmin = ylo95, xmax = yhi95, color = variable, shape = variable, group = variable)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  ggstance::geom_linerangeh(position = ggstance::position_dodgev(.25)) +
  ggstance::geom_pointrangeh(aes(xmin = ylo90, xmax = yhi90), size = 1.25, position = ggstance::position_dodgev(.25)) +
  ylab("") +
  scale_color_manual(values = cbPalette[-3]) +
  xlab("Effect of Border Proximity on Perceived Cultural, Security,\nand Environmental Outcomes (Std. Deviations)") +
  theme_bw() + 
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"), 
        legend.title = element_blank())

ggsave(here("figures", "Fig7.pdf"), height = 5, width = 7)



# Regression Tables -------------------------------------------------------
Texreg(list(mod_protectTrad, mod_calmOrder, mod_enviroProtect),
  custom.model.names = c("Protect Traditions", "Law \\& Order", "Prosperity over Environment"),
  custom.gof.rows = list("Weights" = c("Yes", "Yes", "Yes"),
                         "Controls" = c("Yes", "Yes", "Yes")),
  caption = "Effects on Outcomes from VOX. Model 1: Effects on Protecting Traditions. Model 2: Effects on Emphasizing Law and Order. Model 3: Effects on Favoring Economic Prosperity Over Protecting the Environment. Sample: Border Region, Travel Minutes $<=$ 30. Data: The Post-Popular Vote Survey (VOX). Regressions include municipality and survey wave fixed effects. Controls include income, gender, age, education level, employment status, and marital status.",
  label = "tab:vox_results",
  file = here("tables", "vox_results.tex")
  )

